home *** CD-ROM | disk | FTP | other *** search
/ Trading on the Edge / Trading On The Edge - CD-ROM Toolkit (Wayzata Technology)(2031)(1994).bin / pc / mac_file / vendor_d / neuralwa / nw2v50 / n2rand.c < prev    next >
Text File  |  1993-08-23  |  12KB  |  347 lines
  1. /* $Id$ */
  2. /* $Log$ */
  3. /* n2rand.c - NeuralWare Standard Random Number Generator */
  4. /* Version 1.1 of 6 March 1990 */
  5. /************************************************************************
  6.  * Copyright(C) 1987-1992 NeuralWare Inc                                *
  7.  * Penn Center West, IV-227, Pittsburgh, PA 15276                       *
  8.  * Telephone: (412) 787-8222    FAX: (412) 787-8220                     *
  9.  *                                                                      *
  10.  * All rights reserved.  No part of this program may be reproduced,     *
  11.  * stored in a retrieval system, or transmitted, in any form or by any  *
  12.  * means, electronic, mechanical, photocopying, recording or otherwise  *
  13.  * without the prior written permission of the copyright owner,         *
  14.  * NeuralWare, Inc.                                                     *
  15.  *                                                                      *
  16.  *                          PROPRIETARY NOTICE                          *
  17.  *                                                                      *
  18.  * This document is the property of NeuralWare, Inc. and contains       *
  19.  * trade-secrets and other proprietary information.  The information    *
  20.  * herein is reserved as proprietary to NeuralWare, and is not to be    *
  21.  * published, reproduced, copied, disclosed, used, or reverse           *
  22.  * engineered without the express written consent of a duly authorized  *
  23.  * representative of NeuralWare.                                        *
  24.  ************************************************************************
  25.  */
  26.  
  27.  
  28. /************************************************************************
  29.  *                                                                      *
  30.  *  Random Number Generator for NeuralWare Products                     *
  31.  *                                                                      *
  32.  ************************************************************************
  33.  05-Dec-90 cck
  34.  Modified n2rand.c and n2rand.h for stand-alone operation and
  35.  easy integration with flash code.
  36.  Moved RRand() and GRand() from n2fl_sum.c to n2rand.c
  37.  */
  38.  
  39. /* Its purpose is to provide a standard random number generator across all  */
  40. /* products that NeuralWare sells and supports.                             */
  41.  
  42. /* This code derived from Knuth, The Art of Computer Programming, Vol2,     */
  43. /* Seminumerical Algorithms; Sedgewick, Algorithms; and Press, Flannery,    */
  44. /* Teukolsky, and Vetterling, Numerical Recipes in C                        */
  45.  
  46. /* The concept is to use one type of generator, a linear congruential       */
  47. /* generator, to fill an array with random numbers. We then shuffle that    */
  48. /* array, and then take the difference between two of those random numbers  */
  49. /* in the array. That value is returned to the user, and also used to       */
  50. /* replace one of the numbers in the random number array.                   */
  51.  
  52. /* NRAND Returns a uniform random deviate between 0 and MAXRAND-1.          */
  53. /* Call NRAND as so: (SL) new_deviate = nrand();                            */
  54.  
  55. /* SRAND Sets seed of random sequence to another part of the sequence.      */
  56. /* Call SRAND as so: SRAND( (SL) mi_random_seed );                          */
  57.  
  58. #include <math.h>
  59.  
  60. #include "n2rand.h" /* defines MAXRAND & others       */
  61.  
  62. static long    rand_array[56];           /* Lookup table of random numbers */
  63. static long    first_gen, second_gen;    /* Result for the LCG             */
  64. static int     table_index, loop;        /* Index into random number table */
  65. static int     first_time_flag = 1;      /* First time thru flag -
  66.                                             initially true */
  67. static int     cursor_1, cursor_2, i;    /* Cursors to chose random #'s from
  68.                                             lookup table of random #'s    */
  69. /*   */
  70. /************************************************************************
  71.  *                                                                      *
  72.  *  set_rand() - Random number seed                                     *
  73.  *                                                                      *
  74.  ************************************************************************
  75.  */
  76.  
  77. #ifdef PROTOTYPING
  78. void set_rand( long seed )
  79. #else
  80. void set_rand( seed )
  81. long seed;
  82. #endif
  83. {
  84.    /* Reset first_time_flag on the first call */
  85.    if ( first_time_flag ) {
  86.       /* Now reset first time flag */
  87.       first_time_flag = 0;                /* now false */
  88.    }
  89.    
  90.    /* Now test & correct seed value */
  91.    first_gen = MSEED - (seed < 0 ? -seed : seed);
  92.    first_gen %= MAXRAND;
  93.    rand_array[55] = first_gen;
  94.    second_gen = 1;
  95.    
  96.    /* Now initialize the shuffle table */
  97.    /* We use a linear congruential generator to set up initial table */
  98.    for ( table_index = 1; table_index <= 54; table_index++)
  99.       {
  100.          i = (21*table_index) % 55;
  101.          rand_array[i] = second_gen;
  102.          second_gen = first_gen - second_gen;
  103.          if ( second_gen < MZ ) second_gen += MAXRAND;
  104.          first_gen = rand_array[i];
  105.       } /* End of for table_index = 1 to 54 */
  106.    
  107.    /* Now randomize the table */
  108.    for (loop=1; loop<=4; loop++)
  109.       {
  110.          for ( table_index=1; table_index <= 55; table_index++)
  111.             {
  112.                rand_array[table_index] -= rand_array[1 +
  113.                                                      (table_index+30) % 55];
  114.                if ( rand_array[table_index] < MZ ) rand_array[table_index]
  115.                   += MAXRAND;
  116.             }
  117.       }
  118.    
  119.    /* Finish up any remaining bits */
  120.    
  121.    cursor_1 = 0;                       /* init cursors into random # table */
  122.    cursor_2 = 31;
  123.    
  124.    
  125. }   /* End of initialization/setup code (if first_time true or reinit) */
  126. /*   */
  127. /************************************************************************
  128.  *                                                                      *
  129.  *  nrand() - NeuralWare Random Number Generator                        *
  130.  *                                                                      *
  131.  ************************************************************************
  132.  */
  133.  
  134. long nrand()
  135. {
  136.    /* Retrieve random number & generate a new replacement value */
  137.    /* Initialize on the first call, if not already initialized by set_rand */
  138.    if ( first_time_flag ) {
  139.       set_rand(0);
  140.       /* Now reset first time flag */
  141.       first_time_flag = 0;                /* now false */
  142.    }
  143.    
  144.    if (++cursor_1 == 56) cursor_1 = 1;     /* rollover so 1 to 55 */
  145.    if (++cursor_2 == 56) cursor_2 = 1;     /* rollover so 1 to 55 */
  146.    
  147.    /* now get diff of 2 random #'s */
  148.    first_gen = rand_array[cursor_1]-rand_array[cursor_2];
  149.    
  150.    /* correct for rollover */
  151.    if ( first_gen < MZ ) first_gen += MAXRAND;
  152.    
  153.    /* replace old random # */
  154.    rand_array[cursor_1] = first_gen;
  155.    
  156.    /* and return random # to user */
  157.    return first_gen;
  158. }
  159. /*   */
  160. /************************************************************************
  161.  *                                                                      *
  162.  *  RRand(), GRand() - real random number generators                    *
  163.  *                                                                      *
  164.  *  RRand - Uniform Distribution                                        *
  165.  *  GRand - Gaussian Distribution                                       *
  166.  *                                                                      *
  167.  ************************************************************************
  168.  */
  169.  
  170. double RRand( lv, hv )
  171. double lv, hv;    /* low / high limits */
  172. {
  173.    double rv;    /* return value */
  174.    
  175.    rv = nrand() * RANDINV;   /* rand() / whatever value is needed */
  176.    
  177.    rv = (hv - lv) * rv + lv;   /* range adjusted */
  178.    
  179.    return( rv );
  180. }
  181.  
  182. double GRand( lv, hv )
  183. double lv, hv;    /* low / high limits */
  184. {
  185.    static int    k=0;      /* k=0 need to caluclate g0 & g1,
  186.                             * k=1 g1 already calculated */
  187.    static double g0,g1;    /* Gaussian(0,1) random variables */
  188.    double u0,u1;     /* uniforms on -1 to 1 */
  189.    double r2;        /* radius from (0,0) to (u0,u1) */
  190.    double x,y;       /* intermediate calculations */
  191.    
  192.    if( k ) {
  193.       k = 0;
  194.       return( 0.5 * (hv - lv) * g1 + 0.5 * (hv + lv) );
  195.    } else {
  196.       k = 1;
  197.       /* find two U(-1,1)'s that are independent and lie within unit circle
  198.        * using NeuralWare's stanard uniform generator
  199.        */
  200.       u1 = 2.0 * (nrand() * RANDINV) - 1.0;
  201.       y = u1 * u1;
  202.       for(r2=2.0; r2>1.0; ) {
  203.          u0 = u1;
  204.          r2 = y;
  205.          u1 = 2.0 * (nrand() * RANDINV) - 1.0;
  206.          y = u1 * u1;
  207.          r2 += y;
  208.       }
  209.       /* now get two Gaussian(0,1) */
  210.       x = sqrt( ( -2.0 * log(r2) )/ r2 );
  211.       g0 = u0 * x;
  212.       g1 = u1 * x;
  213.       return( 0.5 * (hv - lv) * g0 + 0.5 * (hv + lv) );
  214.    }
  215. }
  216. /*   */
  217. #ifdef TEST_RAN
  218. /************************************************************************
  219.  *
  220.  *  Test Routine for Random Number Generator
  221.  *
  222.  ************************************************************************
  223.  */
  224.  
  225. /* The following code is a simple test harness for the random number        */
  226. /* generator. This was derived from Sedgewick, Algorithms, Addison-Wesley,  */
  227. /* 1984.                                                                    */
  228.  
  229. #include "usermath.h"
  230.  
  231. #define R_MAX 2000           /* This implies up to 2000 values in histogram */
  232.  
  233. SREAL chi_square( N, r )
  234. UL    N, r;
  235.  
  236. /*
  237.   This routine calculates the sum of the squares of the frequencies
  238.   of occurence of each random index produced via nrand. The sum is
  239.   then scaled by the expected frequency of occurence, and normalized
  240.   by the size of the sequence.
  241.   */
  242.  
  243. /*
  244.   N is the total number of random numbers we wish to generate
  245.   r is the range we wish to produce random number in; The numbers
  246.   generated range from 0 thru r-1.
  247.   
  248.   A good fit is indicated by how close this 'Chi-Square' is to the
  249.   expected value, i.e. for 10000 random numbers ranging between 1 &
  250.   100, a 'good' Chi-square would be around 100 in value
  251.   */
  252.  
  253. {
  254.    static long f[R_MAX];
  255.    static long i, index;
  256.    static SREAL temp, t1, t2;
  257.    
  258.    /* First clear out the histogram array */
  259.    for (index=0; index<R_MAX; index++) f[index] = 0;
  260.    
  261.    /* Now accumulate a histogram of times a given index is produced */
  262.    for (i=0; i<N; i++) {
  263.       
  264. #if 1                      /* set this 1 to test Neuralware generator */
  265.       
  266.       t1 = nrand();
  267.       t2 = r;
  268.       index = (t1*t2)*RANDINV;
  269. #endif
  270.       
  271. #if 0                      /* set this 1 to test alternate generator */
  272.       t1 = random();
  273.       t2 = r;
  274.       index = (t1*t2)/2147483647.0;
  275. #endif
  276.       
  277.       if ( index >= R_MAX || index < 0 ) {
  278.          printf("ERROR - Chi-square, index = %d\n",i);
  279.          index = 0;
  280.       };
  281.       f[index] = f[index] + 1;
  282.    }
  283.    
  284.    /* Now calculate the Chi Square of those frequencies */
  285.    /* We use floating point calculations to see any patterns in low order bits */
  286.    t1 = 0.0;
  287.    t2 = r;
  288.    
  289.    for (i=0; i<=r-1 ; i++)
  290.       t1 = t1 + (f[i]-(N/r))*(f[i]-(N/r));
  291.    
  292.    temp = ((t2*t1)/N);
  293.    
  294.    return temp;
  295.    
  296. }
  297.  
  298. #ifdef __ZTC__                            /* if Zortech C */
  299. #include <time.h>
  300. #endif
  301.  
  302. main()
  303. #define SAMPLE_SIZE 1000
  304. #define RANGE       100
  305. #define NTRIALS     10
  306. {
  307.    NINT     i;            /* loop index */
  308.    SREAL temp[NTRIALS];   /* temporary store */
  309.    SL    seed;            /* seed for random number generator */
  310.    
  311. #ifdef __ZTC__                            /* if Zortech C */
  312.    time_t start, finish;  /* data structures to measure run time of rand funct */
  313. #endif
  314.    
  315.    /* First reseed random number generator */
  316.    seed = 0;
  317.    set_rand(seed);
  318.    
  319.    /* Second, start the loop for Ntrials of random numbers */
  320.    
  321. #ifdef __ZTC__                            /* if Zortech C */
  322.    time (&start);
  323.    for (i=0; i< NTRIALS; i++) {
  324.       temp[i] = chi_square(SAMPLE_SIZE, RANGE);
  325.    }
  326.    time(&finish);
  327. #else
  328.    for (i=0; i< NTRIALS; i++) {
  329.       temp[i] = chi_square(SAMPLE_SIZE, RANGE);
  330.    }
  331. #endif
  332.    
  333.    /* Third, print out the duration per call & CHI square values */
  334.    
  335. #ifdef __ZTC__                            /* if Zortech C */
  336.    printf("time per Random call %2.10f seconds\n",
  337.           (difftime(finish, start)/SAMPLE_SIZE));
  338. #endif
  339.    
  340.    for (i=0; i< NTRIALS; i++)
  341.       printf("Trial # %d    CHI Square = %f\n",i,temp[i]);
  342.    
  343.    
  344.    
  345. }
  346. #endif /* ifdef TEST_RAN */
  347.